About Issues Tutorials Documentation </span>

Tutorial 3: Intro to Quantum Chemistry

This tutorial shows how to select a quantum chemistry method, visualize orbitals, and analyze the electronic wavefunction.



In [ ]:

    
%matplotlib inline
import numpy as np
from matplotlib.pylab import *
try: import seaborn  #optional, makes plots look nicer
except ImportError: pass

import moldesign as mdt
from moldesign import units as u

I. Build and minimize minimal basis set hydrogen

A. Build the molecule

This cell builds H₂ by creating the two atoms, and explicitly setting their positions.

Try editing this cell to:

Create HeH⁺
Create H₃⁺
Change the atoms' initial positions



In [ ]:

    
atom1 = mdt.Atom('H')
atom2 = mdt.Atom('H')
atom1.bond_to(atom2,1)
atom2.x = 2.0 * u.angstrom

h2 = mdt.Molecule([atom1,atom2], name='H2', charge=0)
h2.draw(height=300)

B. Run a hartree-fock calculation

The next cell adds the RHF energy model to our molecule, then triggers a calculation.

Try editing this cell to:

Change the atomic basis
Get a list of other available energy models (type mdt.models. and then hit the [tab] key)



In [ ]:

    
h2.set_energy_model(mdt.models.RHF, basis='3-21g')
h2.calculate()

print('Calculated properties:', h2.properties.keys())

print('Potential energy:', h2.potential_energy)

C. Visualize the orbitals

After running the calculation, we have enough information to visualize the molecular orbitals.



In [ ]:

    
h2.draw_orbitals()

D. Minimize the energy

Here, we'll run a quick energy minimization then visualize how the hydrogen nuclei AND the atomic wavefunctions changed.



In [ ]:

    
minimization = h2.minimize(frame_interval=1, nsteps=10)
minimization.draw_orbitals()

II. Analyzing the wavefunction

The wavefunction created during QM calculations will be stored as an easy-to-analyze python object:



In [ ]:

    
wfn = h2.wfn
wfn

A. Molecular orbital data

First, let's examine the molecular orbitals. The overlaps, fock matrix, coefficents, and density matrix are all available as 2D numpy arrays (with units where applicable).

We'll specifically look at the "canonical" orbitals that result from Hartree-Fock calculations.



In [ ]:

    
mos = wfn.orbitals.canonical
mos

MOs are, of course, a linear combination of AOs:

\begin{equation} \left| \text{MO}_i \right \rangle = \sum_j c_{ij} \left| \text{AO}_j \right\rangle \end{equation}

The coefficient $c_{ij}$ is stored at mos.coeffs[i,j]



In [ ]:

    
mos.coeffs

Most MO sets are orthogonal; their overlaps will often be the identity matrix (plus some small numerical noise)



In [ ]:

    
mos.overlaps

By definition, the fock matrix should be orthogonal as well; the orbital energies are on its diagonal.



In [ ]:

    
matshow(mos.fock.value_in(u.eV), cmap=cm.seismic)
colorbar(label='fock element/eV')
title('Fock matrix')

The MolecularOrbitals class also offers several methods to transform operators into different bases. For instance, the overlap method creates an overlap matrix between the AOs and MOs, where olap[i,j] is the overlap between MO i and AO j: \begin{equation} \text{olap[i,j]} = \left\langle MO_i \middle| AO_j \right \rangle \end{equation}



In [ ]:

    
olap = mos.overlap(wfn.aobasis)
olap

Various other matrices are available from this this object, such as the two-electron Hamiltonian terms:



In [ ]:

    
matshow(mos.h2e.value_in(u.eV), cmap=cm.inferno)
colorbar(label='2-electron hamiltonian term / eV')
title('2-electron hamiltonian')

B. Individual orbitals

You can work with inidividual orbitals as well. For instance, to get a list (in order) of our four atomic orbitals (i.e., the basis functions):



In [ ]:

    
aos = wfn.orbitals.atomic
aos[:]

Let's grab the lowest orbital and examine some of its properties:



In [ ]:

    
orb = aos[0]
print('Name:', orb.name)
print('Energy:', orb.energy)

Orbital objects also give you access to various matrix elements:



In [ ]:

    
ha_1s = aos[0]
hb_1s = aos[3]

print('Overlap between 1s orbitals: ', ha_1s.overlap(hb_1s))
print('Fock element between 1s orbitals', ha_1s.fock_element(hb_1s))

C. Basis functions

An object representing the wavefunction's basis functions is available at wfn.aobasis



In [ ]:

    
wfn.aobasis

It stores a list of AtomicBasisFunction objects:



In [ ]:

    
basis_function = wfn.aobasis[0]
print('Name:', basis_function.name)
print('Angular quantum number:', basis_function.l)
print('Azimuthal quantum number:', basis_function.m)
print('Centered at:', basis_function.center)

Each basis function is defined as a linear combination of "primitive" 3D Gaussian functions:



In [ ]:

    
basis_function.primitives

And these primitives can themselves be examined:



In [ ]:

    
primitive = basis_function.primitives[0]
print(primitive)
print("Coeff:", primitive.coeff)
print("Alpha:", primitive.alpha)



In [ ]: